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SUMMARY 

Results of a calculation of an optimized truncated scarfed nozzle 
were compared. The truncated scarfed nozzle was designed for an exit 
Mach number of 6.0, l.e., the Mach number at the last nozzle 
characteristic Is 6.0, with an external flow Mach number of 5.0. The 
nozzle was designed by the Rao method for optimum thrust nozzles 
modified for two-dimensional flow and truncated scarfed nozzle 
applications. This design was analyzed using a shock-fitting method 
for two-dimensional supersonic flows. 

Excellent agreement was achieved between the design and analysis. 
Truncation of the lower nozzle wall (cowl) revealed that there Is an 
optimum length for truncating the cowl without degrading the nozzle 
performance. Truncation of the nozzle cowl past this optimal length 
should be analyzed In trade-off studies for thrust loss versus gross 
vehicle weight. 

Plots of the oblique shock wave equations were also Identified 
which will allow computation of slip line angle, dynamic pressure 
coefficient, or ambient Mach number for various specific heat ratios. 
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NOMENCLATURE 


a speed of sound 

A coefficient In finite difference equations and cross-sectional area 
AR area ratio, A e /At 
C characteristic curve 

C thrust coefficient - 

2 pV 2 A 

f function 
I Integral 
L length, Inches 
m mass flow rate 
M Mach number 
P static pressure 
P total pressure 
q dynamic pressure 

Q coefficient In finite difference equations 
r radial distance 

R coefficient In finite difference equations 
Rg gas constant 
s distance 

S coefficient In finite difference equations 
t time 

T temperature and coefficient In finite difference equations 
thrust 

u velocity In x-dlrectlon 
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V 


velocity In y-dlrectton 
x Cartesian coordinate 

y Cartesian coordinate 

z axial distance 
Greek Symbols 

Y ratio of specific heats 

6 slip line angle, first variational operator, and coefficient In 
two-dimensional flow equations 

e oblique shock wave angle 

X Lagranglan multiplier and slope of characteristic curve 

V Mach angle 

v Prandtl-Meyer expansion angle 
<|> control surface angle 
p density 

0 oblique shock wave angle 
Subscripts 
a ambient 
e exit 

PM Prandtl-Meyer 
s shock 
si slip line 
t throat 
T total 

+ left-running characteristic 
right-running characteristic 
o streamline 


vl 



INTRODUCTION 


In the design of jet aircraft, It Is desirable to design for 
optimum thrust while restricting the length of the exhaust nozzle for 
overall vehicle weight reduction. In recent years, Interest In 
hypersonic and supersonic cruise vehicles has grown and presently 
studies are being performed to evaluate the feasibility of such 
vehicles. Exhaust nozzles designed for vehicles such as the high-speed 
civil transport (Orient Express) will have to be optimized for thrust 
while conforming to weight limits and other constraints Imposed by 
particular missions. 

Early attempts to design optimum thrust nozzles focused on the 
truncation of perfect nozzles. A perfect nozzle was designed for a 
uniform exit velocity. This nozzle was then truncated at a given area 
ratio. This process was repeated for designs with different uniform 
exit velocities. Each truncation then provided a nozzle with a 
different length for the desired area ratio. Examination of the 
thrust versus length or weight considerations provided the basis of 
selection of the desired optimal design. 

Later attempts to design optimal nozzles Involved the analysis of 
a family of assumed shapes. A shape such as a parabolic contour was 
assumed with a prescribed Initial expansion angle and exit lip angle. 
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The angles were varied to provide a family of nozzle shapes. These 
nozzles were then analyzed to determine the performance and the optimal 
shape was therefore determined. 

A method for designing an optimum thrust nozzle was developed by 
Rao (1958) In which variational calculus and the method of 
characteristics are used to compute the contour. The Rao optimization 
technique Is unique In that It does not depend on a preselected family 
of mathematical functions for the contour and proceeds to find the 
optimal wall contour within that family. The Rao method determines the 
unique contour which yields the maximum thrust from the Infinite array 
of possible contours. It does not necessarily determine a contour 
which can be specified by a simple mathematical equation. 

The Rao method was Incorporated Into a computer program by 
Nickerson (1982). The code computes a nozzle contour of optimum thrust 
(1) for a given length, or (2) for a given area ratio, or (3) for both 
fixed length and exit plane radius. The Rao method was originally 
formulated for axlsymmetric nozzles with circular throats as 
Illustrated In Fig. 1. However, because of the renewed Interest In 
high-speed aircraft the method was modified to compute two-dimensional, 
flat rectangular throat nozzles which can be used In ramjet engine 
applications. Such a nozzle is shown In Fig. 2. The method was also 
applied to designs where the lower nozzle wall (cowl) Is terminated at 
the point where the last characteristic that emanates from the upper 
nozzle wall (ramp) intersects the cowl. This defines a two- 
dimensional, nonsymmetrlc nozzle as shown In Fig. 2. The flat throat 
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INTERNAL 

FLOW 



FIG. 2. - SCHEMATIC OF A TWO-DIMENSIONAL RAO NOZZLE WITH FLAT-THROAT 
MOODIFICATION 
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FIG. 3. - SCHEMATIC OF A TWO-DIMENSIONAL TRUNCATED RAO SCARFED NOZZLE 
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modification assumes an Initial expansion of uniform flow followed by 
a Prandtl-Meyer expansion at the beginning of the turning section. 

Significantly reduced vehicle weights can be obtained with two- 
dimensional, nonsymmetrlc nozzles and further weight reduction can be 
achieved by truncating the nozzle cowl upstream of the last nozzle ramp 
characteristic. However, because of the truncation, the design of the 
nozzle contour will be affected by the external flow as illustrated In 
Fig. 3. If the nozzle flow Is underexpanded, an oblique shock wave 
and slip line will form at the exit because of the Interaction of the 
external and Internal flows. 

The incorporation of external flow considerations Into the 
optimization of exhaust nozzle contours defines the problem that Is 
studied In this thesis. The objective of this work Is to present a 
design method that may be used to compute optimized truncated scarfed 
nozzles and to present results of a parametric Investigation that 
employs this method. 



ANALYSIS 


Rao Method 

A procedure for computing axl symmetric, optimum thrust exhaust 
nozzle contours was developed by Rao (1958). This method defines a 
nozzle contour for an Ideal gas with constant specific heats In an 
Isentroplc flow. For this method, a control surface Is defined at the 
exit of the nozzle. The thrust Is maximized such that the flow and 
length of the nozzle are fixed constraints. The solution of the 
external problem and the flow within the nozzle are then found by 
using the method of characteristics. This procedure was recently 
converted Into a computer program by Nickerson (1982). 

Figure 4 is a schematic of the top portion of a two-dimensional, 
axlsymmetrlc nozzle contour with the characteristic net and control 
surface used for this analysis superimposed upon It. Line CE, In 
Fig. 4 describes the control surface having an angle of Inclination 
<t> to the axis. Figure 5 Is a differential element of the control 
surface showing flow across It. Calculating the mass flow through 
this differential element yields (1) 

dm - pv s1n($ - 0) dA. (1) 

where dA ■ 2ir r ds and ds = dr/sin <J». To obtain the mass flow 
crossing the control surface, the differential mass flow rate In 
Eq. (1) Is Integrated along line CE 
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FIG. 4. - SCHEMATIC OF A RAO NOZZLE WITH CHARACTERISTIC NET AND 
CONTROL SURFACE 
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» - f c p* 2 ” r dr - 

The linear momentum flux In the z-dlrectlon Is: 

momentum flux - pv 2 9 r dr. (3) 

Computing the thrust on the nozzle yields 

*•- l[ [<P - Pa> * P - 2 SM 4 ITH r 

The length of the nozzle Is computed from 

L ■ z + f cot 4> dr. (5) 

c j c 

For a fixed throat contour the length of the nozzle to point C, z c . Is 
also fixed by the application of the entrance flow character and the 
solution of the flow equations. Hence the length constraint becomes 

f cot <|> dr - constant. (6) 

J C 

Utilizing the Lagranglan multiplier method developed In Hildebrand 
(1976) and reviewed In Appendix A, the problem can now be reduced to 
maximizing the following Integral: 


where 


I 



+ x 2 f 2 + Xjf^) dr 



>a > 


+ pv 


slnU - 9) 1 
sin <t> J r 


(7) 


( 8 ) 
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* .. s1n<4> - 0) „ 

f 2 ‘ P tf ~ sTm r <9> 

f 3 - cot (10) 

Variational calculus principles reviewed In appendix B can now be 
applied to Eq. (7). Maximizing I requires that the first variation 
of the Integral must be equal to zero. In turn, this will lead to an 
expression for the control surface and the flow conditions along It. 

The variation of I depends upon the optimization required. If the 
thrust Is maximized such that the length is constant but the area 
ratio varies, then the upper limit of the Integral Is a variable and 
contributes to the variation of I. The variation of I does not 
depend on the lower limit. This can be Illustrated by considering 
that the Integral between limits of points C and E can be expressed as 
the sum of two Integrals. One Integral for the region between points 
C and D and the other Integral for the region between points D and E. 
Since the throat contour Is fixed, and hence, so Is the resultant flow 
the variation of the Integral with the limits of points C and E 
vanishes. Examining all possible variations of the quantities In I, 
the following set of possibilities are derived: 

(1) SC, SM, 50, and S<J> are zero along CD, since this portion of 
the flowfleld Is fixed by the specification of the throat contour. 

(2) SM, S9, and S4> are nonzero along DE. 

(3) At E, Sze Is zero, but Srg Is nonzero. 

(4) Since M and 6 are continuous along CE, then SD, although 
nonzero, does not enter Into the first variation of I. 



Performing the first variation of the Integral I, results In 


pf* r 


61 * 0 


f l M + X 2 f 2 M + X 3 f 3 M J 6M + ^ f l e + X 2 f 2 0 + X 3 f 3 e; ' S0 


,) 


+ [ f 1 + ^ 2^2 + ^3^3 | S4> I dr + Sr^ 1^2 * ^2^2 ^2^3 


♦ 


■9 


9> 


L 


( 11 ) 


where the subscripts M, 0, and 9 denote partial derivatives and 6 
Is the variational operator. Since the variations In 6H, 60, 6$, and 
6r|: are arbitrary, each of their coefficients must be equal to zero, 
thus, 

(12) 


f . + X-f, + X,f, ■ 0 

'm c m j 


f l Q + X 2 f 2 Q + x 3 f 3 0 = 0 

f i + x 2 f 2 + x 3 f 3 - 0 

9 9 9 


(13) 

(14) 


along line DE, and 


fl + X 2 f 2 + X 3 f 3 - 0 at E. 


(15) 


If radial and axial distances at point E are both specified, the upper 

limit Is fixed and the last term In Eq. (11) disappears from the 

variational solution. Since f, = f, - 0, these terms may be 

J M J 0 

eliminated from Eqs. (12) and (13) and the two resulting equations are 
combined to give 
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The density, pressure, and velocity can be considered functions of Mach 
number for Isentroplc flow and can be computed from Isentroplc flow 
relations e.g., see Shapiro (1953) or John (1969). These relations are 


p T 



(17) 



(18) 


v 




(19) 


Note, the static parameters have been nondlmenslonal Ized with respect 
to the total conditions. Utilizing Eqs. (17) to (19), It follows that 


2 


_ pv 

p M ■ - H 

(20) 

P H - -pM 

(21) 

V 

/ \ 

(22) 

m! 1 + M 2 ) . 



(Equations 20 and 21 are nondlmenslonal Ized to total conditions and 

P M = d(p/p T )/dM, p M = dp/dM, and v M = dv/dM). The partial derivatives 

f, and f, are determined from Eqs. (8) and (9) to give: 

'M 


'M 


2 2 2 
pv 2pv s1n(<t> - 6) cos 6 v Mp s1n(<t> - 0) cos 0 

M u f v _ i A " sin 4> 

M 


^1 + M 2 js1n <t> 


(23) 
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and 



s 1 n<4> - 9) 
sin $ 


M 


£V 




(24) 


with p and p nondl mens Iona 1 Ized to total conditions. Slmlllarly, 

the partial derivatives f, and f- are found to be: 

'© ^0 


'6 


-r pv 2 [s1n(<|> - 9) sin 0 + cos(<|> - 0) cos 0 

sin <j> 


(25) 



2 

pv cos(<t> - 9) 
sin <(> 


(26) 


Substituting Eqs. (23) to (26) Into Eq. (16) and rearranging yields: 


£ 

M 


+ 


2p s1n(4> - 9) cos 9 
M fl + *-=-!■ M 2 ! sin <|) 


Mp s 1 n(4> - 9) cos 9 
sin 


COS (<J) - 0) 


s i n (4> - 6) 
sin 0 




X 



0) sin 0 + 


COS(<J) 


0) cos 0 


or 



15 


-p f 1 + Y ~ 1 M 2 j cos <<() - 0) sin <j> + 2p sin(<|> - 0) cos 0 cos(<t> - 0) 


p sin (<|> - 0) + p s1n($ - 0) cos($ - 0) cos 0 


Mp sin 2 (<|> - 0) sin 0 M ( 1 + *-^-1 M 2 


Combining Eqs. (17) and (18) results In the following expression for 
density 


p* - p (i + 3L i J - . 


This Is then substituted Into the relation above and the result Is 
rearranged to give: 

2 1 


M 


- I cos (<(> - 0) - s1n($ - 0) cos 0 cos($ - 0) 

sin* 1 ($ - 0) sin 0 L 


+ sin 2 (<|> - 0) sin ©j 


Utilizing the following trigonometric Identities In this relation 

s1n(<x±3) * sin a cos (3±cos a sin & 

COS (a±|3) = COS a cos |3±s1n a sin |5 


(27) 

(28) 


and solving yields 


M 


1 


sin (<|> - 0) 

This can be expressed In terms of the Mach angle as follows 

p * s 1 n(<(> - 0) * sin y, 

and hence 

= 0 + y. 


(29) 


( 30 ) 
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This shows that the control surface Is a left-running characteristic 
and Eq. (30) Implies that 


^ - tan(9 + y) along DE. 


By Eq. (13) 


2 f. 


Introducing Eq. (30) into Eqs. (25) and (26) and applying the 


trigonometric identities gives: 


r p v cos(9 - u) 
sin $ 


r r pv cos 
2 0 sin <|) 


Substituting these relations Into Eq. (32) and solving for \ 2 gives 

X . -v . cQ ii e ^ K ) . (3: 

2 COS y 

Equation (14) can now be solved for \ 3 

- f l + x 2 f 2 


The partial derivatives f, , f ? and f, can readily be written as: 

♦ ♦ <t> 

f __ pv 2 r cos 6s1n_e (35) 

0 sin 2 4» 

f. = pvr (36) 

<t> sin <j> 
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f, = -=V • (37) 

J <j> sin 

Substituting Eqs. (35) to (37) and (33) Into (34) yields: 

X- = pv 2 r cos 0 sin e - - - c - pvr sin 0 
3 v cos p v 

or 

* 3 = pv 2 r S i„ 9 [ c°s 6 COi_M_^COi(9 - y ) j (38) 

Employing the trigonometric Identity brings 

X 3 = -pv 2 r s1n 2 e tan p. (38) 

Equations (30), (33), and (38) now yield the appropriate 
conditions to calculate the control surface. The constants X 2 and 
X 3 can be evaluated from the properties at point D. According to 
Eqs. (30) the control surface must coincide with a left-running 
characteristic. If this equation Is satisfied, then It follows that 
the compatibility equation for a left-running characteristic must also 
be satisfied along the control surface, line DE. The compatibility 
equation for axl symmetric, Isentroplc, irrotatlonal flow as found In 
Shapiro (1953) Is 


d 6 - 


'M " 1 rIM Slf1 U Slf1 0 Hr fi 

— dM + - .iI/a : \ dr * 0 


■^M 2 


r s 1 n (0 + p) 


( 39 ) 


It can be shown that Eq. (39) Is satisfied by Eqs. (33) and (38). 
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In order for the nozzle to be optimum for a fixed length 

constraint, the condition In Eq. (15) must be satisfied at point E. 

Utilizing Eqs. (8) to (10), (33), and (38), Eq. (15) becomes: 

IY_ _ s .2 sin u cos 0 J v cos(0 + u) sin u 
L (p - V + pV sin(© ♦ V ) r J p pv since . J> r 


2 2 

- rpv sin 6 tan p cot(0 + p) * 0, 


which Is the same as 


P - Pa 
1 2 

2 pv tan p 


-2 cos u cos 0 2 cos(0 + u ) - 

s1n(0 + p) + s1n(0 + p) _ “ c 


2 

sin 0 cot(0 + p). 


The right-hand side of this equation may now be reduced by use of 
trigonometric identities. The final result Is 


P - p 

sin 20 * -j 2 — > at point E. (40) 

2 pv* tan p 

Equations (33) and (38) for X2 and X3 are used In conjunction 
with the characteristic equations to define the flowfleld and which 
relate the variables r, z, M, and 0 are: 



dr = tan(0 + p) dz 

along a left-running characteristic and 


(41) 


(42) 


d0 



(43) 


dr = tan(0 - p) dz 


(44) 
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FIG. 6. - CHARACTERISTIC PLOT FOR TWO-DIMENSIONAL SUPERSONIC 


IRROTATIONAL FLOW 
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along a right-running characteristic, where the Mach angle p Is 
defined by: 



The equations developed above are schematically shown In Fig. 6. 

For two-dimensional Cartesian geometries Eq. (38) Is independent 
of r and becomes 

X 3 = -pv 2 sin 2 © tan p. (46) 

Equations (30) and (33) remain the same as for axlsymmetrlc flow and 
along with Eq. (46) provide two equations In the two unknowns M and 
6 and hence yield a constant M and © along the last 
characteristic. The character of flow is therefore that of a simple 
wave and this Is In agreement with flow in the turning section of a 
two-dimensional Cartesian geometry nozzle. 

Scarfed Truncation Point Computation 
A schematic of the nomenclature used In the scarfed computer 
subroutine for computing the flow properties at a truncation point on 
the nozzle cowl are shown In Fig. 7. Truncation of the nozzle to 
point F allows the external flow to affect the mechanics of the 
solution. Since the nozzle Is assumed to be underexpanded, l.e., the 
Internal static pressure Is greater than the ambient static pressure, 
a Prandtl-Meyer expansion fan and an oblique shock wave emanate from 
point F. Since the flow properties will differ In regions 2 and V In 


21 


REGION AFFECTED 



FIG. 7. - SCHEMATIC OF THE NOMENCLATURE USED IN SCARFED 
TRUNCATION POINT SUBROUTINE 
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Fig. 7, a slip line also forms at point F from the expansion of the 
Internal flow. These flow phenomena are Incorporated Into the optimum 
nozzle program. 

Figure 8 Is a computer flow diagram of the scarfed subroutine. 
This subroutine computes the phenomena which occur at the cowl 
truncation point. Data Is Input for the ambient and nozzle Mach 
numbers M a and M] respectively, specific heat ratios and static 
pressures. A check Is performed on the ambient Mach number to 
determine whether Its value Is greater than 1.0. If M a Is 
supersonic, the static pressure behind the Prandtl-Meyer expansion Is 
Initialized by setting It equal to the ambient static pressure. If 
M a Is subsonic, the static pressure behind the Prandtl-Meyer 
expansion Is Initialized by computing the average of the ambient and 
nozzle static pressures. The flow regions computed by the scarfed 
subroutine are numbered as shown In Fig. 7. 

Total pressures for the ambient flow and flow regions 1 and 2 
along with the Prandtl-Meyer expansion angle v, at state 1 are now 
computed from the isentroplc flow equations developed In standard 
compressible flow texts, e.g. see John (1969) or Shapiro (1953). These 
relations are: 





(47) 


and 




(M 2 - 1) - tan -1 



- 1 ' 


V ■ 


(48) 



FIG. 



8. - COMPUTER FLOW DIAGRAM FOR THE SCARFED TRUNCA- 


TION POINT SUBROUTINE 
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The total pressure at state 2 is set equal to the total pressure at 
state 2', because the total pressure across the slip line Is constant. 
The Mach number at state 2 is next computed from Eq. (47) along with 
the corresponding Prandtl-Meyer expansion angle. The total Prandtl- 
Meyer expansion angle is then computed by calculating the absolute 
value of the difference in expansion angles at states 1 and 2. 

The next step in the scarfed subroutine is to perform a test on 
the ambient Mach number to establish Its value. For ambient Mach 
numbers less than 1.0, the subroutine returns to the main Rao program. 
For ambient Mach numbers greater than 1.0, the slip line angle, 6, and 
oblique shock wave angle, 6, are computed. The slip line angle Is 
computed from 


COt 6 


tan 6 


<Y - 1)M* 

M? s1n 2 0 - 1 


- 1 


(49) 


and the oblique shock wave angle from 

P, 2yM? Sln 2 e - ( Y - 1) 

ST ■ V ,1 ' <5 ° 

Both expressions are developed In the Ames Tables (1953). 

The absolute value of the difference of the total Prandtl-Meyer 
expansion angle and the oblique shock wave angle, phi, Is computed and 
checked for convergence. A difference of 0.00001 between successive 
passes was used as the convergence criteria In this study. If the 
convergence test fails, the subroutine computes the derivatives of 
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Eqs. (47) to (50) and these are used In Newton’s method for Iterative 
solution of simultaneous equations as developed In Press et al . (1986) 
Newton's method consists of computing the difference of the Initial 
guesses and the derivatives of functions which define the guesses. 


The derivatives In the matrix above are defined as: 


1 In a 

2-by-: 

dv 

dS 

dM 

de 

dp 2 

dp 2 

dM^ 

W 

: above 

are < 


dp 2 

dMl 


L Y ♦ 1 


M,p 


2 K 2 


1 + 2-y-i M* 


dv 

dM 


-1 




! ( ! + M* 


dS 

sin & cos & 

1 

de " 

Pa 

(E - 1) 


1 


yM^ - 5 


- 1 


- [2 Y M^ - (y - 1) - (y + 1)?] 


-1 


- :L f 1 [<Y + ♦ <Y - 1>] 


-1 


(Note In the above £ = p 2 /p a -> 
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2 

dp, 4y M^p a sin 0 cos 0 

c aad 

d0 Y, + 1 

d 

This square matrix can be solved by the determinant method. The value 
of the determinant, phlc, Is based upon the values of the flow 
quantities computed for each Iteration. The difference In total 
Prandtl -Meyer expansion angle and oblique shock wave angle, phi Is now 
divided by phlc and this value Is then used to correct the Initial 
guess of the state 2 static pressure. This Is expressed In equation 
form as follows 


'n-1 


Phi 

Phlc 


n-l 


1*2 n. 


The computational procedure outlined Is now repeated until convergence 
Is realized. The oblique shock wave angle, state 2 Mach number and 
static pressure are then used In the Rao program for constructing the 
right-running characteristic at the truncation point. The computer 
listing of this subroutine Is contained In Appendix C. 


Slip Line and Oblique Shock Wave Computation 
The oblique shock wave and slip line that occur at the cowl 
truncation point may be curved or straight dependent upon the pressure 
difference between the ambient and nozzle flowflelds. In general, the 
oblique shock wave and slip line will be curved because of the 
additional expansion which takes place In the Internal nozzle 
flowfield. If the oblique shock wave is curved, the strength of the 
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shock will vary from point to point. The flow downstream of the shock 
Is only Isentroplc along streamlines and Is rotational because of 
gradients In the entropy and stagnation properties normal to the 
streamlines. For Isentroplc flows, the entropy Is constant along 
streamlines. 

Zucrow and Hoffman (1977) develop the governing equations for 
steady adiabatic flow of an Invlscld compressible fluid, which are 
used to derive the characteristic equations for a rotational 
Isentroplc flow. The continuity, momentum, and speed of sound (used 
In place of energy), equations are: 


V • (pv) - 0 

p 0T + 7p * 0 


Dp 2 Dfi_ 
Dt 3 Dt 


0 . 


Expressing Eqs. (51) to (53) In Cartesian coordinates for two- 
dimensional flow yields: 

pu x + pv y + up x + vp y + ^ " 0 
puu x + pvu y + P x - 0 

puv x + P vv y + p y = 0 
up x + vp y - a 2 up x - a 2 v Py - 0 

where 6*0 for planar flow, and 6 = 1 for axl symmetric flow. 


(51) 

(52) 

(53) 


(54) 

(55) 

(56) 

(57) 
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The characteristic and compatibility equations which result from 
Eqs. (54) to (57) for a two-dimensional isentropic rotational 
supersonic flow developed in Zucrow and Hoffman (1977) are: 

(a) Characteristic Equations 

= \ = 77 (streamline) (58) 

o u 



(HxJ ± = x ± " tan(0± l i) 

(b) Compatibility Equations 


V7^2 

pv 2 


dp ± d6 + 8 


(Mach lines) 


(59) 


pv dv + dp = 0 

(along streamlines) 

(60) 

2 

dp - a dp = 0 


(along streamlines) 

(61) 

sin 0 dx 




± 

yM cos(0±p) 

= 0 

(along Mach lines). 

(62) 


Figure 9 shows a plot of the characteristics for steady two-dimensional 
supersonic rotational flow. The characteristic and compatibility 
equations derived for the Mach lines remain the same as those equations 
derived for irrotational flow. However, two additional equations are 
obtained. These are the streamlines (counted twice) and the 
compatibility relations along the streamlines. The compatibility 
equations are the Bernoulli equation and the speed of sound relation. 

The finite difference equations derived by Zucrow and Hoffman for 
two-dimensional steady isentropic rotational flow equations are: 

Ay„ = V 

■'o o o 

R Av + Ap = 0 
oo o 


( 64 ) 




X 

FIG. 9. - CHARACTERISTIC PLOT FOR STEADY TWO-DIMENSIONAL SUPER- 
SONIC ROTATIONAL FLOW 
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Ap - A Ap = 0 
r 0 0 r 0 

(65) 

Ay ± * \ AX ± 

(66) 

Q ± Ap ± A0 ± + S ± Ax ± - 0 

(67) 


(68) 

pv 

(69) 

20 

O 

II 

T> 

< 

(70) 

r. 8 sin 0 

* " yM cos(0±p) 

(71) 


where C_ denotes a right-running characteristic, C + denotes a 
left-running characteristic and C 0 denotes a streamline as shown In 
Fig. 9. 

The numerical algorithms for applying Eqs. (63) to (71) are based 
upon a modified Euler predictor-corrector method. In the predictor 
step, each of the coefficients In the finite difference equations are 
calculated at the known Initial points. In the corrector step, average 
values of the four primary dependent variables v, 0, p, and p are 
computed along each three of the characteristics, and the coefficients 
of the finite difference are calculated by employing those average 
properties. This numerical Integration algorithm Is termed an average 
property method. 

The steady two-dimensional supersonic rotational flow equations 
derived above have to be Incorporated into the Rao computer program to 
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complete the construction of the flowfleld between the external oblique 
shock wave and the slip line. Figure 10 provides an Indication of the 
changes Incorporated Into the Rao computer program. The top portion, 
Fig. 10(a), shows left and right-running characteristics emanating from 
the nozzle cowl. The bottom portion, Fig. 10(b) Is a schematic of the 
mesh used to calculate the rotational flow. The truncation point on 
the nozzle cowl corresponds to grid point 11 on the mesh. The flow 
conditions at this grid point were calculated by the scarfed routine 
discussed previously. 

The flow properties at grid point 21 are computed by assuming that 
the flow properties at this point are defined by right-running flow 
characteristics from the Internal flow (which Is discussed below) with 
an assumed Mach number. The program then Iterates using the continuity 
relation between the left-running characteristic from point A to 
point 11 and the right-running characteristic from point A to point 21 
to establish the final value of the conditions at point 21. Once the 
flow properties are determined at grid point 21, the properties at 
grid points 11 and 21 are used to find the flow properties at grid 
point 23 using the oblique shock wave relations and the method of 
characteristics for rotational isentroplc flow. The conditions at 
points 21 and 23 are then Iterated upon until convergence is realized 
In the pressure loss across the oblique shock wave. The flow 
properties at grid point 22 are computed by taking an average of the 
properties at grid points 21 and 23. The flow properties at following 
mesh points are subsequently computed. 
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FIG. 11. - FINITE DIFFERENCE NETWORK FOR TWO-DIMENSIONAL 
SUPERSONIC ROTATIONAL FLOW 
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Interior Point Computation 

As the calculation proceeds, a numerical algorithm Is required 
for calculation of an Interior point. The finite difference network 
used In this analysis Is Illustrated In Fig. 11 . Figure 12 Is a 
schematic of the unit process for determining an Interior point. 

Points 1 and 2 are Initial-data points on the C_ and C + Mach lines 
and line 12 connects these points. The rearward running streamline 
C 0 Intersects line 12 at point 3 and the properties at point 3 are 
computed by linear Interpolation. 

The following four equations are obtained by writing Eqs. (63) 
and ( 66 ) In finite difference form In terms of points 1 , 2 , 3 , and 4 , 
and from the equation of the diagonal line 12. 

y 3 ' X o X 3 = y 4 - X o X 4 <72) 

>4 - \. x 4 * y 2 " X + x 2 <73) 

y H - X_x 4 = y 1 - X_x 1 (74) 

y 3 ■ x 1 2 x 3 s — ^ 12 X 2 (75) 

where X] 2 , X + , X_, and Xq, are the slopes of line 12, the left- 
running characterl Stic, the right-running character! stlc, and the 
streamline, respectively. Equations (72) to (75) are solved 
simultaneously for X 3 , y 3 , X 4 , and y 4 . The quantities v, e, p, 
and p are then calculated by linear Interpolation as discussed 
above. The Euler predictors for an Interior point are: 


X + = tan(0 2 + p 2 ) 


(76) 
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FIG. 12. - SCHEMATIC OF THE UNIT PROCESS FOR AN INTERIOR POINT 
COMPUTATION 
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X_ * tan(e 1 -| ij). 


(77) 


The value of Xq depends upon the flow angle at point 3, which 
Is not known until point 3 Is located and the corresponding flow 
properties have been calculated by linear Interpolation. An Iterative 
procedure must be employed for determining the location of point 3 
during each pass through the modified Euler predictor-corrector 
algorithm. The Initial estimate for 03 Is obtained for the Euler 
predictor as follows: 


Equations (76) and (77) are then solved simultaneously for X 3 and 
y 3 , and the corresponding ©3 Is found by linear Interpolation. The 
procedure may then be repeated using the new value ©3 to obtain 
Improved values of X3, y 3, and 63. The procedure Is repeated until 
the change In the values of X 3 and y 3 diminish below the specified 
tolerance. 

The compatibility equation, Eq. (67) which Is valid along Mach 
lines, may be written In finite difference form as: 

Q + P 4 + e 4 - T + (79) 

Q_P 4 + e 4 - T_ , (80) 

where 


T + * - S^(x 4 - X 2 ) + Q + P£ + (81) 

T_ = -S_(x 4 - Xj) + Q_p 1 - © 1 . (82) 


For the Euler predictor the following equations are derived: 
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» 2 V 2 

VTT 

pi v f 


5 sin 6„ 


s + * y 2 M 2 cos (© 2 + ji 2 > 


(83) 


6 sin 6 


y^ cos(© 1 + jij ) 


(84) 


These equations are now solved simultaneously for P4 and ©4. 

Writing Eqs. (64) and (65), the compatibility equations valid 
along streamlines. In finite difference form gives 

R o v 4 + P4 " R o v 3 + P3 " T 0l 
P4 " A q p 4 " p 3 " A o p 3 * T 0 o 


(85) 

( 86 ) 


For the Euler predictor 

R o = P 3 V 3 • A o * a 3 <87> 

where a * a(p, p). Equations (85) and (86) are now solved for V4 
and P4. This completes the application of the Euler predictor 
algorithm. The computational equations derived for steady two- 
dimensional Isentroplc rotational supersonic flow are summarized as 
follows: 

y-5 “ X « X 1 m y A ~ \>> X * 


r 3 V*3 J 4 V4 

y- - X x. - y 0 - X x« 
J 4 + 4 J 2 + 2 

- V 4 “ - X _ x i 


(73) 

(74) 


y 3 ” X 12 x 3 " y 2 " X 12 x 2 


(75) 
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Q p. + 0- - T 
Q_P 4 - e 4 * T _ 
R 0 W 4 * P« • T o, 

p 4 ' A o p 4 * T o 2 


(79) 

(80) 

(85) 

( 86 ) 


Improved values for the location of and the corresponding flow 
properties at point 4 may be computed by repeating the steps outlined 
above for the Euler predictor algorithm, but employing average values 
of v, 0, p, and p along each characteristic. Hence, 

X + « tan<0 + + p + ) (88) 

where 


P + 


P2 * P4 

2 


0 

+ 


9 2 * 6 4 
2 


v 

+ 


V 2 + V 4 
2 


P 2 + P4 

Pj. * 5 


(89) 


For a simple thermodynamic system 

d » 

and It follows 


a(p, p) 


a 


+ 


a(p + , p + ) 


Then, 


(90) 

(91) 



a - sin 
+ 



(92) 


The coefficients X_, Xq, Q + , Q_, S + , S_, R 0 , and A 0 are determined 
In a similar manner. 

The speed of sound, a, now has to be determined from Eq. (91). 
Assuming a perfect gas, Eq. (91) takes the form 
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a 


2 


, YB 


P 


(93) 


Correspondingly, the temperature T, Is given by 

T = T(p, p) 


and for a perfect gas 



(94) 


(95) 


The pressure p and density p are used to calculate the temperature 
T, Mach number M, and the speed of sound a, In a separate subroutine. 
By calculating the fluid properties In this manner, the different unit 
processes can be determined by the method of characteristics, which are 
Independent of the equations of state for the fluid. 

Iteration and convergence are controlled on the axis of symmetry. 
The ratio sin ©2/y2 In the coefficient S Is approximated by the 
ratio sin ©i /y i for the predictor. For the corrector, that ratio Is 
based upon average values of 6 and y, which are nonzero. 


External Shock Wave Point Computation 
Figure 13 Illustrates schematically the unit process for 
determining an external shock wave point. Properties are known at 
point 1 on the downstream side of the shock wave, and point 2 Is a 
known point on the previous right-running Mach line 12. The free- 
stream flow properties Moo, v®, p*,, and poo are also known. Point 4, 
on the downstream side of the shock wave, Is the external shock wave 
point. This point Is located at the intersection of the shock wave 
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FIG. 13. - SCHEMATIC OF THE UNIT PROCESS FOR AN EXTERNAL SHOCK 


WAVE POINT COMPUTATION 
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from point 1 and the left-running Mach line 24. The slope of shock 
wave 14 Is computed from the average value of the two wave angles e ] 
and E 2 at the points 1 and 4, respectively. These wave angles are 
functions of the free-stream Mach number M a and the flow turning 
angles 0 ] and 64 . 

To determine the location of, and the properties at point 4, an 
Iterative procedure Is employed. A value of the wave angle 64 Is 
assumed at point 4. The finite difference equation for the shock 
wave 14 Is 

1 + e 4 > \ 

y 4 - y , - V* 4 - V - tan I — 2 — I <x 4 _ V <96) 

and the equation of the left-running Mach line 24 Is 

y 4 - y 2 - V x 4 ' x 2 > - (97) 

The location of point 4 Is computed by solving Eqs. (96) and (97) 
simultaneously for X 4 and y 4 - Equation (79), the single 
compatibility equation valid for left-running Mach lines, must also be 
satisfied for this computation. 

Incorporation of Rotational Flow 

The equations developed for the Interior point and external shock 
wave point computation have to be Incorporated Into the Rao program. 

The flow diagram In Fig. 14 Illustrates how this can be accomplished. 

A point on the slip line Is computed as previously described. A 
comparison Is now performed to determine If the ambient static pressure 
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FLOW DIAGRAM FOR ROTATIONAL 
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Is less than the nozzle static pressure. If the ambient static 
pressure Is greater than the nozzle static pressure, the subroutine 
returns to the main program (this case would depict a shock wave 
forming In the Interior of the nozzle and this analysis will fail for 
this case). If the ambient static pressure Is less than the nozzle 
static pressure, the computation continues. 

Next, the Interior point Is examined to determine whether It Is 
the first point after the scarfed truncation point. If It Is the first 
point, the routine computes a shock point using the external shock wave 
point routine. A shock loss Is then computed and compared to the 
pressure behind the shock wave. If this loss agrees with the pressure 
behind the shock wave additional Interior points are considered. These 
points are moved and readjusted Into arrays. The routine then advances 
to the next characteristic and the computation returns to the main 
program. 

On the other hand. If the point computed on the slip line Is not 
the first point after the scarfed truncation point, an Interior point 
computation Is performed. If this Interior point Is upstream of the 
first slip line point (see Fig. 10), this point Is discarded and the 
next point Is used. If this point Is downstream of the scarfed 
truncation point, all Interior points are calculated. A shock point Is 
subsequently computed and the computation continues as described above. 
The computation of the kernel continues downstream until the contour Is 
defined. 
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Seagull Analysis 

The design of a truncated scarfed nozzle can be analyzed using a 
shock fitting method for two-dimensional supersonic flows. This shock 
fitting method referred to as Seagull, was developed by Salas (1976). 
Seagull can be used to analyze complex two-dimensional or axl symmetric 
supersonic Invlscld flows of a perfect gas. Seagull Is limited to 
analyzing flows whose axial component of Mach number remains supersonic 
(hyperbolic equations of motion are used In this analysis and for 
subsonic flows these equations become elliptic). This code can be used 
to compute flowflelds produced by a single duct or several ducts that 
are merged. Jets and plumes can also be analyzed by Seagull. This 
code will be used In this thesis for a comparison with the results from 
the Rao method. 


Flow Correlations 

Manipulation of Eqs. (49) and (50) produces a functional 
relationship for the external and nozzle Mach numbers, oblique shock 
wave and slip line angles, ambient and nozzle static pressures, and 
dynamic pressure coefficient. 

Rearranging Eqs. (49) and (50) yield the following relations: 



(98) 


<p 2 - p 1 )/q 1 
2 - (p 2 - p ] )/q 2 - 


tan 8 


(99) 
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These relations exhibit a mathematical relationship of a linear nature, 
l.e, having the form y = mx + b, and can be graphically presented. 

This can provide a simple graphical representation of the flow 
conditions at the truncated cowl lip point. 



RESULTS AND DISCUSSION 


Flow Correlations 

Equations (99) and (100), the relations relating the slip line and 
oblique shock wave angles, ambient and nozzle static pressures, are 
plotted In Figs. 15 and 16. 

Figure 15 Is a plot of dynamic pressure coefficient times specific 
heat ratio versus oblique shock wave angle for constant values Mach 
number. Straight line curves result for this plot and this Illustrates 
the linear nature of Eq. (99). For a given dynamic pressure 
coefficient, ambient specific heat ratio, and ambient Mach number, the 
oblique shock wave angle can be determined or for any of three 
previously mentioned flow parameters, a fourth can be determined. 

A plot of oblique shock wave angle versus dynamic pressure 
coefficient for constant values of slip line angle Is shown In Fig. 16. 
The curves In this plot were constructed for slip line angles of 1.0°, 
2.0°, 10.0°, 25.0°, and 50.0°. These curves exhibit a hyperbolic 
behavior due to the cotangent trigonometric function In Eq. (100). As 
the slip line angle Is Increased the curves become progressively linear 
and tend towards Infinity as the dynamic pressure coefficient Is 
Increased. 


46 


47 



FIG. 15. - PLOT OF Sin 2 OF OBLIQUE SHOCK WAVE ANGLE VERSUS 
DYNAMIC PRESSURE COEFFICIENT 



100 


80 





DYNAMIC PRESSURE COEFFICIENT, P 2 "P 1 /q 1 


FIG. 16. - PLOT OF DYNAMIC PRESSURE COEFFICIENT VERSUS 
COTANGENT OF OBLIQUE SHOCK HAVE ANGLE 
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Truncated Scarfed Nozzle Computation 

A truncated scarf nozzle was calculated for a design Mach number 
of 6.0, l.e, the Mach number at last nozzle characteristic Is 6.0 with 
an external flow Mach number of number of 5.0, and external and 
Internal specific heat ratios of 1.4. An area ratio of 90.0, nozzle 
ramp length of 300.0 inches, nozzle cowl length of 15.0 Inches, nozzle 
throat Mach number of 1.24, nozzle throat temperature of 2500 °R, 
along with ambient static pressure of 0.150 psla and 100 psla for the 
nozzle throat total pressure complete the Input used for this design. 
Pressures and coordinates computed In this analysis are 
nondlmenslonal Ized with respect to the nozzle throat conditions. 

A nontruncated scarfed nozzle was also designed to compare the 
effect of truncation. The nontruncated nozzle was computed for the 
same design Mach number as the truncated nozzle (the external flow does 
not affect the nozzle computation, refer to Fig. 2). The nontruncated 
nozzle computation requires approximately 10 seconds of cpu time on an 
IBM 3033 mainframe computer, while the truncated scarfed case requires 
approximately 20 seconds of cpu time. 

Figure 17 Is a plot of nozzle cowl length versus nozzle thrust 
coefficient. The scarfed nozzle designed for a Mach number of 6.0 was 
selectively truncated for nozzle cowl lengths ranging from 50 Inches to 
12.5 inches. A nozzle thrust coefficient of 1.6661 Is computed for 
the nontruncated case and Is Invariant with cowl truncation down to 
20.0 inches. Truncating the cowl shorter than 20.0 inches will cause 
the nozzle thrust coefficient to drop sharply as shown In Fig. 17. A 
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nozzle cowl length of 20.0 Inches corresponds to point where the last 
nozzle ramp characteristic Intersects the cowl surface as Illustrated 
In Fig. 2. 

The nontruncated nozzle and truncated scarfed nozzle were then 
analyzed with the Seagull code. These computations were performed on 
the Cray X-MP supercomputer and required approximately 30 seconds of 
cpu time per case. 

Figures 18 to 22 are plots of flow path geometry, static pressure, 
axial velocity, axial Mach number, and Mach number contours for the 
nontruncated nozzle. Figures 19 and 20 display smooth distributions of 
static pressure and axial velocity plots, where the Mach number 
distribution In the plot of Fig. 21 shows a wave In the lower wall Mach 
number. The wave can be attributed to the Prandtl-Meyer expansion fan 
which occurs at the beginning of the turning section past the throat 
section and the Intersection of the two families of characteristics. 
This slight compression wave Is clearly shown In Fig. 22 In the Mach 
number contour plot. 

The analysis of the truncated scarfed nozzle case are presented In 
Figs. 23 to 28. Figure 23 Is a plot of the flow path geometry with the 
cowl truncated at 15.0 Inches. The solid line extending past this 
point Is a boundary drawn by the Seagull program. Figure 24 contains 
the static pressure plot of the truncated scarfed nozzle and shows an 
abrupt drop In the static pressure distribution on the lower wall at 
the truncation point. The lower wall static pressure becomes equal to 
the freestream static pressure slightly beyond the truncation point. 
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FIG. 20. - WALL AXIAL VELOCITY PLOT OF NONTRUNCATED NOZZLE 
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A 1.2400 
B 1.4406 
C 1.6411 
D 1.8417 
E 2.0423 
F 2.2429 
6 2.4434 
H 2.6440 
I 2.8446 



MACH NUMBER 


FIG. 22. - MACH NUMBER CONTOUR PLOT OF NONTRUNCATED NOZZLE 
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FIG. 23. - FLOW PATH GEOMETRY OF TRUNCATED SCARFED NOZZLE 


STATIC PRESSURE 



FIG. 24. - HALL STATIC PRESSURE PLOT OF TRUNCATED SCARFED 


NOZZLE 
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Figure 25 displays the computed axial velocity distribution. As 
exhibited In the static pressure distribution, the axial velocity drops 
to the freestream velocity after the truncation point. The Mach number 
distribution also displays a similar behavior and shows that the 
truncated scarfed nozzle does not Internally expand as much as the 
nontruncated nozzle. This clearly Illustrates the effect of the 
external flow on the expansion of the Internal flow as shown In 
Fig. 26. 

Figures 27 and 28 are plots of the oblique shock wave and slip 
line and Mach number contour for the truncated scarfed nozzle. 

Figure 27 displays the oblique shock wave and slip line emanating from 
the cowl truncation point. The slip line Initially diverges from the 
Internal flow and then converges as the effect of the external flow 
Increases with nozzle length. The Mach number contour plot In Fig. 28 
clearly Illustrates the Interaction of the external and Internal flows. 
The slip line Is depicted by the coalescing of contour lines In the 
center of the plot and the strength of the oblique shock wave Is shown 
In the bottom portion of the figure. Also, the Prandtl-Meyer 
expansion fan Interacts with the slip line In this case Instead of the 
nozzle cowl wall In the nontruncated case. 

Figure 29 Is a comparison of nozzle ramp wall Mach Number 
distributions for the truncated scarfed nozzle. The Rao design and 
Seagull analysis compare very favorably except for the corner region 
where the compression wave occurs. This Is Illustrated by the jump In 
Mach Number at approximately 2.0 Inches as shown In Fig. 29. 
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FIG. 25. - WALL AXIAL VELOCITY PLOT OF TRUNCATED SCARFED 


NOZZLE 
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STATION LOCATION 


FIG. 26. - WALL AXIAL MACH NUMBER PLOT OF TRUNCATED 


SCARFED NOZZLE 
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A 1.2400 
B 1.4465 
C 1.6529 
D 1 .8594 
E 2.0658 
F 2.2723 
G 2.4788 
H 2.6852 



MACH NUMBER 


FIG. 28. - MACH NUMBER CONTOUR PLOT OF TRUNCATED SCARFED NOZZLE 
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FIG. 29. - COMPARISON OF NOZZLE RAMP WALL MACH NUMBER 


DISTRIBUTIONS 


CONCLUDING REMARKS 


Results of a calculation of an optimized truncated scarfed nozzle 
were compared. The scarfed nozzle design showed less expansion than 
the nontruncated case and this was primarily due to the external flow 
affecting the Internal flow expansion within the nozzle. 

The comparison of thrust coefficient versus nozzle cowl length 
revealed that truncation of the cowl will affect the overall 
performance of an exhaust nozzle. This comparison demonstrates that 
there Is an optimal cowl length In which truncation can be performed 
without degrading the overall nozzle performance. Truncation of the 
cowl past this optimal length should be analyzed in trade-off studies 
for thrust loss versus gross vehicle weight. 

Plots of the oblique shock wave equations will allow the 
calculation of oblique shock wave and slip line angles, dynamic 
pressure coefficient, or ambient Mach number for various specific heat 
ratios. This will allow a designer to compute these flow quantities 
directly from the plots presented without use of the compressible flow 
equations. 

The truncated scarfed nozzle method developed In this thesis can 
be extended in the future. The method can be modified to compute an 
optimum truncated scarfed nozzle with an oblique shock wave forming 
Internally and simple chemical kinetics could also be incorporated 
Into this method. 
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APPENDIX A 


LAGRANGIAN MULTIPLIERS 

Situations occur In many engineering applications In which a 
function f Is to be maximized or minimized dependent upon variables 
which are not Independent of each other, but are interrelated by one 
or more constraint conditions. 

The general case may be Illustrated by the problem of maximizing or 
minimizing a function f(x,y,z), so that f(x,y,z) = relative max or 
min subject to two constraints of the form 


g(x,y,z) = 0, 

(A. 1 ) 

h(x,y,z) = 0, 

(A. 2) 


where g and h are not functionally dependent, so that the 
constraints are neither equivalent nor compatible. It Is assumed that 
the functions f, g, and h have first partial derivatives everywhere 
In a region which Includes the desired critical point. 

The function f Is an extremem at P(x 0 , y 0 , z 0 ) If the linear 
terms In Taylor expansion of f about P are zero. This condition 
can be written In the form 

f x dx + fy dy + f 2 dz - 0 at (x 0 , y 0 , z 0 >. (A. 3) 

However, the Increments dx, dy, and dz are not Independent, so It 
cannot be concluded that f x , f y , and f z must vanish separately at 
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P. Since g and h are constant throughout the field, the 
differential of each must be zero at the point P, so that 

g x dx + g y dy + g z dz - 0 at (x 0 , y 0 , z 0 ) (A. 4) 

h x dx + hy dy + h z dz = 0 at (x 0 , y 0 , z 0 ). (A. 5) 

Equations (A. 4) and (A. 5) are linear equations In dx, dy, and dz. 
These equations can be uniquely solved for two of the differentials In 
terms of a third If and only if g and h are functionally 
Independent. By making use of Eqs. (A. 4) and (A. 5) two of the 
differentials (say dx and dy) can be eliminated from Eq. (A. 3). 

This results In an equation of the form F(x,y,z) dz = 0, in which the 
one remaining differential can be arbitrarily assigned. In essence 
the conditions In Eqs. (A. 4) and (A. 5) along with the condition 
F(x,y,z) » 0 constitute a set of three equations in the three unknowns 
x 0 . Vo. and z o- 

An alternative method of frequent usefulness. Is based on the fact 
that one may multiply the equal members of Eqs. (A. 4) and (A. 5) by 
arbitrary constants, say, \] and \2» respectively, and add the 
results to Eq. (A. 3) to obtain the requirement that 
(f x + X]g x + X 2 h x ) ^ + (f y + \29y ^ X 2 hy) dy 

+ (f z + Xi g z + X 2 h z ) dz * 0 (A. 6) 

at (x 0 , y 0 , z 0 ). This expression is valid for all values of X] and 
X 2 * The latter which can now be determined so that the coefficients 
of two of the differentials are zero. If the preceding statement were 
not true. It would follow that 
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2 

- 0 
x 

and hence g and h would be functionally dependent, so that the two 
constraints would be either equivalent or Inconsistent. Assuming 
that \] and X2 have been so determined, then the remaining 
differential can be arbitrarily assigned, so that Its coefficient also 
must vanish. Thus, the three equations 

f x + Vx + X 2 h x * 0 

f y + Vy + * 2 h y - 0 at <x o’ V V (A ' 7) 

f z + X l9z + X 2 h z - 0 
along with the conditions 

g - 0 , h = 0 at <x 0 , y 0 , z 0 ) (A. 8) 

comprise five equations In the five unknown quantities x 0 , y 0 , z 0 , Xj , 
and 

The parameters \] and X2 are referred to as Lagrangian 
multipliers. If these parameters were eliminated from Eq. (A. 6), the 
result would be the relation F * 0 obtained by eliminating dx, dy, 
and dz from Eqs.(A. 3 ) to (A. 5 ). 

The relations in Eq. (A. 7 ) provide the necessary conditions 
£x = 0, £y « 0, * 0 so that the "auxiliary function" 

5 - f + Xig + X2h (A. 9 ) 

attains a relative maximum or minimum at (x 0> y 0 , z 0 ) when no 
constraints are Imposed. 
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The more general case would be If f were to be maximized subject 
to the restraints gi =0, g£ - 0,..., g n - 0. The auxiliary function 
would then be 

K = f + ^191 + ^292 •••■*■ ^n9n» (A. 10) 

where X] * .... X 2 are unknown constants. The necessary conditions 
for rendering £ a relative maximum or minimum with no constraints 
could be written and the solution could be subsequently determined. 

Further explanation of the use of Lagranglan multipliers can be 
found In Hildebrand (1965 and 1976). 



APPENDIX B 


CALCULUS OF VARIATIONS 

The calculus of variations Is a mathematical procedure employed 
to determine one or more functions, subject to certain conditions, so 
as to maximize or minimize a specific definite Integral, whose 
Integrand depends upon the unknown function or functions and/or 
certain combination of their derivatives. The brief discussion of the 
theory of calculus of variations which follows Is presented In greater 
detail In Hildebrand (1976). 

To Illustrate this method an example Is considered. A case Is 
considered In which we wish to maximize or minimize an Integral of the 
form 


I 


I 


F(x,u,u' ) dx. 


(B.l) 


subject to the conditions 

u(a) =• A, u(b) = B, (B.2) 

where a, b, A, and B are specified constants. The function F has 
continuous second-order derivatives with respect to Its three 
arguments. It Is further required that the unknown function u(x) 
possess two derivatives everywhere In the Interval (a, b). 

Functions which have two derivatives In (a, b) and which take on 
the prescribed end values are admissible and out of all admissible 
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functions, a function (or functions) must be selected which makes I a 
maximum. Assuming there Is one function u(x) having the property 
above, a one-parameter family of admissible competing functions which 
Includes u(x), Is first written as 

u(x) + en(x) 

where n<x) Is any arbitrary chosen twlce-dlfferentlable function which 
vanishes at the end points of the Interval (a, b), l.e, 

n<a) - n(b> - 0, (B.3) 

and where e Isa parameter which Is constant for any one function In 
the set but which varies from one function to the next. The Increment 
en(x>, representing the difference between the varied function and the 
actual solution function, Is called a variation of u(x>. 

Replacing u(x) by u(x) + en(x) In I yields another Integral 
that Is denoted by 1(e), 

r b 

1(e) » F(x,u + en.u' + en') dx. (B.4) 

J a 

It then follows that 1(e) takes on a maximum value when e ■ 0, that 
Is, when the variation of u Is zero. Hence, It must follow that 

» 0 when e - 0. (B.5) 

de 

The partial derivatives of F are assumed to be continuous with 
respect to Its three arguments and this Implies that the derivative 
dF/de exists. Hence, 1(e) under the Integral sign may be 
differentiated to obtain 
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dl(e) f raF<x,u + en .u 1 + e n‘ ) _ . 3F(x,u + en.u 1 ± en') 

~Tc J a L 3<u + en> n 3<u’ + en') n J ax - 


(B.6) 


By setting e = 0, an expression for the condition <B.5) Is obtained, 
viz. , 

5I ' <0) ■ J* [m 

a 

In which it Is written that F » F(x,u,u‘) and therefore the partial 
derivatives 3F/3u and 3F/3u ' have been formed with x, u, and u' 
treated as independent variables. 

Transforming the Integral of the second product In Eq. (B.6) via 
an Integration by parts yields 

b Ij pb 

fu 7 dx - [|[jT n<x>] - 


3F 




pU 

d_ 

faF.) 

dx 

J a 

l 3u 'J 


as a consequence of Eq. (B.3). Equation (B.6) now becomes 

n<x) dx - 0 


pU 


(—] 

3F 


dx 

Uu* 1 

3u 

J a 

L 

\ J 



n<x) dx 


(B.7) 


It can be proven that, Eq. (B.7) Is true for any function n(x) 
which Is twice differentiable in (a, b) and zero at the ends of that 
interval. Consequently, the coefficient of n<x) in the integrand 
must be zero everywhere in (a, b), so that the condition 
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d_ f 9F_^ 3£ 

dx l 8u ' J ' 9u 


0 


<B.8) 


must be satisfied. This Is called the Euler equation associated with 
the problem of maximizing or minimizing the Integral In Eq. (B.l) 
subject to Eq. (B.2). 

Recalling that F, and hence Its partial derivatives, may depend 
upon x both directly and indirectly, through the Intermediate 
variables u(x> and u'(x), It Is deduced from the chain rule that 


dM 9M 9M du 9M d 2 u , R 
dx " 9x + 9u dx + 9u' dx 2 

where the function M may be Identified with F or with one of its 
partial derivatives. Thus, In particular, we can employ M - 9F/9u' 
to write the Euler equation in an expanded form as 


2 

c d_U r du /r 

u 1 u 1 dx 2 + r uu' dx + ^ r xu' 


- F u> - 0 


(B. 10) 


The expanded form shows that when F u < u « = 9 2 F/9u 1 2 is zero, the 
equation is In fact a differential equation of second order In u, 
subject to the two boundary conditions u(a> = A and u(b) = B. 
From Eq. (B.8) it follows that 
9F 


9u 


- * constant 


when 


F = F( x , u ' ) , 


(B.ll) 


so that, when F does not Involve u explicitly, the first-order 
equation 9F/9u‘ = constant comprises a "first Integral" of the Euler 
equation. It can also be shown that Eq. (B.8) implies that 
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u' - F * constant when F * F(u,u'), (B.12) 

oU 

so that the first Integral Is also available when F does not Involve 
x explicitly. 

The equation above can be generalized to maximize or minimize an 
Integral subject to a constraint condition. To maximize or minimize 
Eq. (B.l) 

pb 

F(x,u,u’> dx ■ max or min, (B.13) 

J a 

where u(x) Is to satisfy the prescribed end conditions 

u(a) - A, u(b> * B, (B.14) 

as before, but that also a constraint condition Is Imposed In the form 

G(x,u,u' ) dx » K, (B.15) 

where K Is a prescribed constant. For this case, the appropriate 
Euler equation Is found to be the result of replacing F In Eq. (B.8) 
by the auxiliary function 

H * F + XG (B. 16) 

where X Isa Lagranglan multiplier and Is equal to an unknown 
constant. This constant will generally appear In the Euler equation 
and In Its solution. This constant has to be determined along with 
the two constants of Integration In such a way that the three 
conditions of Eqs. (B.15) and (B.16) are satisfied. 



APPENDIX C 


SCARFED TRUNCATION POINT SUBROUTINE 


A computer routine for determining the flow properties at the 
truncation point Is presented of a two-dimensional nozzle Is 
presented. The following equations are used In this analysis 




M 


y/y-i 


(C.l) 




(M - 1) - tan 


- 1 Vm 2 - 


cot 6 = tan e 


(y + 1 )Mj 


2(M 2 s1n 2 0 - 1 


- 1 


(C.2) 


(C.3) 


p 2 2yM 2 s1n 2 e - (y - 1) 


'1 


Y + 1 


(C.4) 


These equations are used to compute static and total pressure, 

Mach number, slip line angle, oblique shock wave angle, and Prandtl- 
Meyer expansion angle at the scarfed truncation point. These flow 
relations can be used to analyze the flow phenomena at the truncation 
point. The truncation point and flow phenomena described previously Is 
schematically shown In Fig. Cl. A computer flow diagram Is presented 
In Fig. C2. The listing of the computer program Is also provided. 
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FIG. Cl. - SCHEMATIC OF TRUNCATED TWO-DIMENSIONAL RAO NOZZLE 


! 
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OOOOIOO c 

0000200 C 
0000300 
0000400 
0000500 
0000600 
0000700 
0000800 
0000900 
0001000 
0001100 
0001200 
0001300 
0001400 
0001500 
0001600 
0001700 
0001800 C 
0001900 C 
0002000 C 
0002100 
0002200 
0002300 
0002400 C 
0002500 C 
0002600 C 
0002700 
0002800 C 
0002900 C 
0003000 C 
0003100 
0003200 
0003300 
0003400 C 
0003500 C 
0003600 C 
0003700 
0003800 C 
0003900 C 
0004000 C 
0004100 
0004200 
0004300 C 
0004400 C 
0004500 C 
0004600 


...THIS ROUTINE COMPUTES THE FLOW QUANTITIES REQUIRED FOR THE SCARFED 
...NOZZLE OPTIMTIZATION 

REAL NU1,NU2<200),M1,M2(200),MA,P1,P2(200),PA,PT1.PT2,PTA 
REAL GM , GMA , DELSL < 200 ) , DELPM (200), THETA ( 200 ) , PHI ( 200 ) , PHIC ( 200 ) 
REAL DNU2DM( 200) , DDLDTH( 200) , DP2DM2 ( 200) ,DP2DTH(200) 

REAL A,B,X,Y,Z,P,G,H 

REAL PR1 , PR2 , PR3 , PR4 , PR5 

READ ( 5, 10, END* 1 ) Ml , MA, PA, GM, GMA 

10 F0RMAT(5(F12 . 5/) ) 

1 1-1 

77 READ (5,11) PI 

11 F0RMAT{F10.5) 

IF ( PI .EQ. 1.0 ) GO TO 30 
IF( MA .GT. 0.0) GO TO 15 
GO TO 20 

...INITIALIZE STATIC PRESSURE AT STATE 2, P2 

15 P2(I)-(Pl+PA)/2. 

GO TO 25 
20 P2(I)-PA 

...CALCULATE THE TOTAL PRESSURE AT STATE 1 

25 PT1-P1*((1.+ (GM-l.)/2. •<M1**2.Q))**(GM/(GM-1.))) 

...CALCULATE THE PRANDTL-MAYER ANGLE NU, AT STATE 1 

NU1-(SQRT ( { GM+1 • ) /*< GM-1 • ) ) ) • ATAN( SQRT( (GM-1 .)/ (GM+1 .) * 

$ (Ml**2 . -1 . ) ) ) - ATAN( SQRT(M1**2.-1.) ) 

PT2-PT1 

...CALCULATE THE MACH NUMBER AT STATE 2 

35 M2(I)»SQRT( ( (PT2/P2(I) ) ••( (GM-1 . )/GM ) -1.) • (2. /(GM-1.)) ) 

...CALCULATE THE PRANDTL-MAYER ANGLE NU, AT STATE 2 

3 NU2 { I ) *SQRT( ( GM+1 .)/( GM-1 .) ) • ATAN( SQRT( ((GM-1.)/(GM 
$+l.))*(M2(I)**2.-l.) ) ) - ATAN( SQRT(M2( I) **2.-1 . ) ) 

...CALCULATE THE TOTAL PRANDTL-MAYER ANGLE, DELPM 


DELPM ( I ) *NU2 ( I ) -NU1 
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0004700 C 
0004800 C 
0004900 C 
0005000 
0005100 C 
0005200 C 
0005300 C 
0005400 
0005500 
0005600 C 
0005700 C 
0005800 C 
0005900 
0006000 
0006100 
0006200 C 
0006300 C 
0006400 C 
0006500 
0006600 C 
0006700 C, 
0006800 C 
0006900 
0007000 C 
0007100 C. 
0007200 C 
0007300 
0007400 
0007500 
0007600 
0007700 
0007800 
0007900 
0008000 
0008100 
0008200 
0008300 
0008400 
0008500 C 
0008600 C. 
0008700 C 
0008800 
0008900 
0009000 
0009100 C 
0009200 C. 


ORIGINAL PAGE !S 

OF POOR QUALITY 


. . .CALCULATE THE AMBIENT TOTAL PRESSURE 

PTA-PA* ( ( 1 .+(GMA-1 . ) /2 . • (MA**2 . ) ) •• (0MA/(0MA-1 . ) )) 

...CALCULATE THE SHOCK NAVE ANGLE, THETA 

THETA(1)-ARSIN( SQRT( ( (P2 ( I ) /PA) • (GMA+1 . )+(GMA-l . ) )/(2. *GMA 
$*MA**2. ) ) ) 

...CALCULATE THE SLIP LINE ANGLE, DELSL 
A" (GMA+1 . ) *MA**2. 

B-2.*; MA**2 . * (SIN(THETA( I) )**2. ) - 1. ) 

DELSL :i)-ATAN( l./( TAN ( THETA ( I ) ) * (A/B - 1.) ) ) 

...CALCULATE PHI FOR ITERATION SCHEME 

PHI(I)-( DELPM(I) - DELSL(I) ) 

...CHECK PHI FOR CONVERGENCE 

IF( ABS( PHI ( I ) ) .LT. 0.000001 ) GO TO 30 

..COMPUTE DERIVATIVES FOR ITERATION SCHEMES 


DNU2DM(I)-SQRT( (H2(I)**2.-1.0) ) /( M2 ( I ) »*3 . • ( (GM-1 . ) /2 . ) 
8+ M2 ( 1 ) ) 

DP2DTH(I)-( 4. *6MA*MA**2. *PA*SIN(THBTA( I) ) •COS(THETA(I) ) )/ 
$ ( GMA + 1.) 

X— 1./( 2 . • SIN ( THETA ( I ) ) •C0S(TH1:TA( I ) ) ) 

Y»( 2 . *MA**2 . *SIN(THETA( I) ) *COS( THETA ( I) ) )/( MA*«2.«SIN( 
$THETA(I) )**2.-l . ) 

2-( 4.*MA»*2.*SIN(THETA(I))*C0S(THETA(I)) )/( 2.+MA**2.*( 

$( GMA+1.) - 2*( SIN ( THETA (I)))** 2.) ) 

OOLDTH ( I ) »SIN ( DELSL ( I ) ) • COS ( DELSL ( I ) ) • ( X+Y+Z ) 

DP2DM2 ( I )--( GM*M2 ( I ) *P2 ( I ) )/( 1 . +( GM-1 . ) »0. B*M2 | I ) ••Z . ) 

..CALCULATE CAP PHI 


PHIC(I)-( DNU2DM( I ) /DP2DM2 ( I ) - DDLDTH ( I ) /DP2DTH( I ) ) 
N-I+l 

ft P2(N)-P2|I) - PHI (I) /PHIC(I ) 

..INDEX COUNTER 
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0009300 C 
0009400 
0009500 
0009600 C 
0009700 C 
0009800 C 
0009900 
0010000 
0010100 
0010200 
0010300 
0010400 
0010500 
0010600 
0010620 
0010700 C 
0010800 C 
0010900 C 
0011000 
0011100 
0011200 
0011300 
0011400 
0011500 
0011600 
0011700 
0011800 
0011900 
0012000 
0012100 
0012200 
0012300 
0012400 
0012500 
EOF 


I-I+l 
GO TO 35 

• • • CONVERT ANGLES TO DEGREES 

30 DELSL ( I ) -DELSL (I) * (1/0.01745239) 

DELPM ( I ) -DELPM ( I ) • ( 1/0 . 01745239 ) 

NU1«NU1*< 1/0.01745239) 

NU2 ( I ) -NU2 { I ) • J 1/0 . 01745239 ) 

PR1-PA/PTA 

PR2-P1/PA 

PR3-PTA/PT1 

PR4-PT1/P1 

PR5-P2(I)/PA 

...WRITE FLOW QUANTITIES COMPUTED 

WRITE ( 6 , 50) NU1 , Ml , PT1 , MA,PTA, I 
50 FORMAT ( 6F10 .5,13,//) 

WRITE (6,60) P 2 ( I ) , M2 ( I ) , DELSL ( I ) , DELPM ( I ) , NU2 ( I ) 

60 FORMAT (5F10. 5 ,/) 

WRITE ( 7 f 70) P2( I) , M2 ( I ) , DELSL ( I ) » DELPM ( I ) * NU2 ( I ) , 

$NU1 ,M1 ,P1 ,PT1 ,MA,PA,PTA,GM,GMA,I 
70 FORMAT { / / # 1 P?« 1 , F10 . 5 , 5X, 'M2- 1 , F10 . 5 , 5X, 'DELSL-* ,F10.5, //, 

$ 1 DELPM- 1 , F9 . 5 , 3X , 'NU2-' ,F10.A # 4X, •NUl« t ,F10.5,//, 'Ml- 1 , 
$F10.5,5X, 'PI-' ,F10.5,5X, *PT1« ' , E10. 5 , // , 'MA-* , F10.5.5X, 'PA-', - 
$F10.5,5X, 1 PTA— • # E10 • 5 , // , •GM-' ,F10.5,5X, 1 GMA-* , F10.5,4X, 'I-MS) 
WRITE ( 8 , 80) PR2 # PR3 ,NU1 # NU2 ( I ) # DELSL ( I ) , MA, PA,M1 , PT1 ,PTA 
80 FORMAT (10F1 1.3) 

WRITE ( 9 , 89 ) PR2 , PR5 , DELSL ( I ) 

89 FORMAT (3F1 1.3) 

99 STOP 
END 
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